home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpkdrv.zip / SQR.FOR < prev    next >
Text File  |  1984-01-06  |  20KB  |  665 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SQRTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SQRTS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER
  15. C
  16. C     TESTS
  17. C        SQRDC,SQRSL
  18. C
  19. C     LINPACK. THIS VERSION DATED 08/14/78 .
  20. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  21. C
  22. C     SUBROUTINES AND FUNCTIONS
  23. C
  24. C     LINPACK SQRDC,SQRSL
  25. C     EXTERNAL SQRBM,SQRFM,SQRRX,SXGEN,SMACH
  26. C     FORTRAN AMAX1,ABS,FLOAT
  27. C
  28. C     INTERNAL VARIABLES
  29. C
  30.       INTEGER LUNIT
  31.       INTEGER CASE,LDX,N,P,PD2,PD2P1,NCASE,JPVT(25)
  32.       INTEGER I,INFO,J,JJ,JJJ,L
  33.       REAL BESTAT,RMAX,RSTAT,XBMAX,XBSTAT,ERRLVL,BSTAT,FSTAT
  34.       REAL BETA(25),QRAUX(25),QY(25),RSD(25),S(25),WORK(25)
  35.       REAL X(25,25),XBETA(25),XX(25,25),Y(25),Y1(25),Y2(25)
  36.       REAL SMACH,TINY,HUGE
  37.       LOGICAL NOTWRT
  38.       LDX = 25
  39.       NCASE = 9
  40.       ERRLVL = 100.0E0
  41.       NOTWRT = .TRUE.
  42.       TINY = SMACH(2)
  43.       HUGE = SMACH(3)
  44.       WRITE (LUNIT,600)
  45.       DO 550 CASE = 1, NCASE
  46.          WRITE (LUNIT,560) CASE
  47.          XBSTAT = 0.0E0
  48.          BESTAT = 0.0E0
  49.          GO TO (10, 100, 170, 240, 300, 330, 380, 430, 470), CASE
  50. C
  51. C        WELL CONDITIONED LEAST SQUARES PROBLEM.
  52. C
  53.    10    CONTINUE
  54.             WRITE (LUNIT,640)
  55.             N = 10
  56.             P = 4
  57. C
  58. C           GENERATE A MATRIX X WITH SINGULAR VALUES 1. AND .5.
  59. C
  60.             PD2 = P/2
  61.             PD2P1 = PD2 + 1
  62.             DO 20 L = 1, PD2
  63.                S(L) = 1.0E0
  64.    20       CONTINUE
  65.             DO 30 L = PD2P1, P
  66.                S(L) = 0.5E0
  67.    30       CONTINUE
  68.             CALL SXGEN(X,LDX,N,P,S)
  69.             IF (NOTWRT) GO TO 40
  70.                WRITE (LUNIT,610)
  71.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  72.    40       CONTINUE
  73. C           THE OBSERVATION VECTOR IS Y = Y1 + Y2, WHERE Y1 IS
  74. C           THE VECTOR OF ROW SUMS OF X AND Y2 IS ORTHOGONAL TO
  75. C           THE COLUMNS OF X.
  76. C
  77.             DO 60 I = 1, N
  78.                Y1(I) = 0.0E0
  79.                Y2(I) = 2.0E0*TINY
  80.                IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)*TINY
  81.                DO 50 J = 1, P
  82.                   X(I,J) = X(I,J)*TINY
  83.                   Y1(I) = Y1(I) + X(I,J)
  84.                   XX(I,J) = X(I,J)
  85.    50          CONTINUE
  86.                Y(I) = Y1(I) + Y2(I)
  87.    60       CONTINUE
  88. C
  89. C           REDUCE THE MATRIX.
  90. C
  91.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  92.             IF (NOTWRT) GO TO 70
  93.                WRITE (LUNIT,610)
  94.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  95.                WRITE (LUNIT,620)
  96.                CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  97.    70       CONTINUE
  98. C
  99. C           COMPUTE THE STATISTICS FOR THE FOWARD AND BACK
  100. C           MULTIPLICATIONS.
  101. C
  102.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  103.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  104. C
  105. C           SOLVE THE LEAST SQUARES PROBLEM.
  106. C
  107.             CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
  108.      *                 INFO)
  109. C
  110. C           COMPUTE THE LEAST SQUARES TEST STATISTICS.
  111. C
  112.             RSTAT = 0.0E0
  113.             RMAX = 0.0E0
  114.             XBSTAT = 0.0E0
  115.             XBMAX = 0.0E0
  116.             DO 80 I = 1, N
  117.                RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
  118.                RMAX = AMAX1(RMAX,ABS(Y2(I)))
  119.                XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  120.                XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
  121.    80       CONTINUE
  122.             RSTAT = RSTAT/(SMACH(1)*RMAX)
  123.             XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
  124.             BESTAT = 0.0E0
  125.             DO 90 J = 1, P
  126.                BESTAT = AMAX1(BESTAT,ABS(1.0E0-BETA(J)))
  127.    90       CONTINUE
  128.             BESTAT = BESTAT/SMACH(1)
  129.             WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
  130.          GO TO 540
  131. C
  132. C        4 X 10 MATRIX.
  133. C
  134.   100    CONTINUE
  135.             WRITE (LUNIT,650)
  136.             N = 4
  137.             P = 10
  138.             PD2 = P/2
  139.             PD2P1 = PD2 + 1
  140.             DO 110 L = 1, PD2
  141.                S(L) = 1.0E0
  142.   110       CONTINUE
  143.             DO 120 L = PD2P1, P
  144.                S(L) = 0.5E0
  145.   120       CONTINUE
  146.             CALL SXGEN(X,LDX,N,P,S)
  147.             IF (NOTWRT) GO TO 130
  148.                WRITE (LUNIT,610)
  149.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  150.   130       CONTINUE
  151.             DO 150 I = 1, N
  152.                DO 140 J = 1, P
  153.                   XX(I,J) = X(I,J)
  154.   140          CONTINUE
  155.   150       CONTINUE
  156.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  157.             IF (NOTWRT) GO TO 160
  158.                WRITE (LUNIT,610)
  159.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  160.                WRITE (LUNIT,620)
  161.                CALL SARRAY(QRAUX,N,N,1,-4,LUNIT)
  162.   160       CONTINUE
  163.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  164.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  165.             WRITE (LUNIT,570) FSTAT,BSTAT
  166.          GO TO 540
  167. C
  168. C        PIVOTING AND OVERFLOW TEST.  COLUMNS 1, 4, 7 FROZEN.
  169. C
  170.   170    CONTINUE
  171.             WRITE (LUNIT,670)
  172.             N = 10
  173.             P = 3
  174.             S(1) = 1.0E0
  175.             S(2) = 0.5E0
  176.             S(3) = 0.25E0
  177.             CALL SXGEN(X,LDX,N,P,S)
  178.             P = 9
  179.             DO 190 I = 1, N
  180.                DO 180 JJJ = 1, 3
  181.                   J = 4 - JJJ
  182.                   JJ = 3*J
  183.                   X(I,JJ) = HUGE*X(I,J)
  184.                   X(I,JJ-1) = X(I,JJ)/2.0E0
  185.                   X(I,JJ-2) = X(I,JJ-1)/2.0E0
  186.   180          CONTINUE
  187.   190       CONTINUE
  188.             IF (NOTWRT) GO TO 200
  189.                WRITE (LUNIT,610)
  190.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  191.   200       CONTINUE
  192.             DO 220 J = 1, P
  193.                JPVT(J) = 0
  194.                DO 210 I = 1, N
  195.                   XX(I,J) = X(I,J)
  196.   210          CONTINUE
  197.   220       CONTINUE
  198.             JPVT(1) = -1
  199.             JPVT(4) = -1
  200.             JPVT(7) = -1
  201.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  202.             IF (NOTWRT) GO TO 230
  203.                WRITE (LUNIT,610)
  204.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  205.                WRITE (LUNIT,620)
  206.                CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  207.   230       CONTINUE
  208.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  209.             CALL SQRRX(XX,LDX,N,P,JPVT)
  210.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  211.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  212.             WRITE (LUNIT,570) FSTAT,BSTAT
  213.          GO TO 540
  214. C
  215. C        25 BY 25 MATRIX
  216. C
  217.   240    CONTINUE
  218.             WRITE (LUNIT,680)
  219.             N = 25
  220.             P = 25
  221.             DO 250 I = 1, 25
  222.                S(I) = 2.0E0**(1 - I)
  223.   250       CONTINUE
  224.             CALL SXGEN(X,LDX,N,P,S)
  225.             IF (NOTWRT) GO TO 260
  226.                WRITE (LUNIT,610)
  227.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  228.   260       CONTINUE
  229.             DO 280 J = 1, P
  230.                JPVT(J) = 0
  231.                DO 270 I = 1, N
  232.                   XX(I,J) = X(I,J)
  233.   270          CONTINUE
  234.   280       CONTINUE
  235.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  236.             IF (NOTWRT) GO TO 290
  237.                WRITE (LUNIT,610)
  238.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  239.                WRITE (LUNIT,620)
  240.                CALL SARRAY(QRAUX,P,P,1,-10,LUNIT)
  241.   290       CONTINUE
  242.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  243.             CALL SQRRX(XX,LDX,N,P,JPVT)
  244.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  245.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  246.             WRITE (LUNIT,570) FSTAT,BSTAT
  247.          GO TO 540
  248. C
  249. C        MONOELEMENTAL MATRIX.
  250. C
  251.   300    CONTINUE
  252.             WRITE (LUNIT,690)
  253.             N = 1
  254.             P = 1
  255.             X(1,1) = 1.0E0
  256.             IF (NOTWRT) GO TO 310
  257.                WRITE (LUNIT,610)
  258.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  259.   310       CONTINUE
  260.             XX(1,1) = 1.0E0
  261.             JPVT(1) = 0
  262.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  263.             IF (NOTWRT) GO TO 320
  264.                WRITE (LUNIT,610)
  265.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  266.                WRITE (LUNIT,620)
  267.                CALL SARRAY(QRAUX,P,